started: Alexey Larionov, 01Mar2016
last updated: Alexey Larionov, 23Apr2019

Summary

Remove genotypes (set to NA) if:
1) genotypes gq < 20 (~12%)
2) genotypes dp < 10 (~18%, +5% after genotypes)
3) genotypes dp > 1000 (~1%)
4) Alt fraction

is > 0.2 in Hom-Ref
is <0.2 or >0.8 in Het
is < 0.8 in Hom-Alt

Remove variants with
1) call_rate < 0.85
(in each: nfe and wecare after the genotypes filtering)
this step removes ~73% of all variants !!!
2) uniform genotypes accross all samples
(created by the above filtering)

In addition
1) rename gt.mx to gt.mx
2) remove unnecessary matrices (gt_num, meta, gq, ref, alt); keep dp and gt_chr for later use

gq 20 filter is set arbitrary; however, it is consistent with what some othersis do
(e.g. see Carson BMC Bioinformatics. 2014 15:125).

A small number of genotypes (~1%) was covered too high to be true (above 1000 coverage).
These are obvious mistakes, and they have been removed too. Arbitrarily the threshold for
max DP was set to 1000 (appr. 10 fold of average coverage).

gt NA rate before filtering ~5% (~3% in wecare and ~12% in NFE respectively)
gt NA rate after filtering ~4% (~3% in wecare and ~6% in NFE respectively)

Input data: 8,275 vars x 739 cases (541 wecare + 198 nfe)
Output data: 1.886 vars x 739 cases (541 wecare + 198 nfe)

start_section

# Time stamp
Sys.time()

## [1] "2019-04-23 15:00:36 BST"

# Clenan-up
rm(list=ls())
graphics.off()

# Base folder
library(knitr)
base_folder="/Users/alexey/Documents/wecare/ampliseq/v04_ampliseq_nfe/s07_filter_genotypes_and_variants"
opts_knit$set(root.dir = base_folder)
options(stringsAsFactors = F,
        warnPartialMatchArgs = T, 
        warnPartialMatchAttr = T, 
        warnPartialMatchDollar = T)
#setwd(base_folder)

# Thresholds for genotypes
min.gq <- 20
min.dp <- 10
max.dp <- 1000
min.call.rate <- 0.85

# add parameters to filter by Alt fraction here?

load_data

load(paste(base_folder, "r01_filter_variants.RData", sep="/"))

gt.mx <- gt_add.mx
rm(gt_add.mx, gt_num.mx, meta.df) # keep gt_chr.mx for TsTv ratio check later

check_data

# List objects
ls()

##  [1] "alt.mx"        "base_folder"   "dp.mx"         "fixed.df"      "gq.mx"         "gt_chr.mx"     "gt.mx"         "max.dp"        "min.call.rate" "min.dp"        "min.gq"        "ref.mx"

# Check sizes
dim(gt.mx)

## [1] 8275  739

dim(ref.mx)

## [1] 8275  739

dim(alt.mx)

## [1] 8275  739

dim(gq.mx)

## [1] 8275  739

dim(dp.mx)

## [1] 8275  739

dim(fixed.df)

## [1] 8275  110

colnames(fixed.df)

##   [1] "CHROM"                   "POS"                     "ID"                      "REF"                     "ALT"                     "QUAL"                    "FILTER"                  "init_AC"                 "init_AF"                 "init_AN"                 "AS_FilterStatus"         "AS_VQSLOD"               "AS_culprit"              "BaseQRankSum"            "ClippingRankSum"         "DB"                      "DP"                      "DS"                      "END"                     "ExcessHet"               "FS"                      "HaplotypeScore"          "InbreedingCoeff"         "MLEAC"                   "MLEAF"                   "MQ"                      "MQRankSum"               "NDA"                     "NEGATIVE_TRAIN_SITE"     "POSITIVE_TRAIN_SITE"     "QD"                      "ReadPosRankSum"          "SOR"                     "SplitVarID"              "VQSLOD"                  "culprit"                 "exac_non_TCGA.AC"        "exac_non_TCGA.AC_AFR"   
##  [39] "exac_non_TCGA.AC_AMR"    "exac_non_TCGA.AC_Adj"    "exac_non_TCGA.AC_EAS"    "exac_non_TCGA.AC_FEMALE" "exac_non_TCGA.AC_FIN"    "exac_non_TCGA.AC_Hemi"   "exac_non_TCGA.AC_Het"    "exac_non_TCGA.AC_Hom"    "exac_non_TCGA.AC_MALE"   "exac_non_TCGA.AC_NFE"    "exac_non_TCGA.AC_OTH"    "exac_non_TCGA.AC_SAS"    "exac_non_TCGA.AF"        "exac_non_TCGA.AN"        "exac_non_TCGA.AN_AFR"    "exac_non_TCGA.AN_AMR"    "exac_non_TCGA.AN_Adj"    "exac_non_TCGA.AN_EAS"    "exac_non_TCGA.AN_FEMALE" "exac_non_TCGA.AN_FIN"    "exac_non_TCGA.AN_MALE"   "exac_non_TCGA.AN_NFE"    "exac_non_TCGA.AN_OTH"    "exac_non_TCGA.AN_SAS"    "exac_non_TCGA.Hemi_AFR"  "exac_non_TCGA.Hemi_AMR"  "exac_non_TCGA.Hemi_EAS"  "exac_non_TCGA.Hemi_FIN"  "exac_non_TCGA.Hemi_NFE"  "exac_non_TCGA.Hemi_OTH"  "exac_non_TCGA.Hemi_SAS"  "exac_non_TCGA.Het_AFR"   "exac_non_TCGA.Het_AMR"   "exac_non_TCGA.Het_EAS"   "exac_non_TCGA.Het_FIN"   "exac_non_TCGA.Het_NFE"   "exac_non_TCGA.Het_OTH"   "exac_non_TCGA.Het_SAS"  
##  [77] "exac_non_TCGA.Hom_AFR"   "exac_non_TCGA.Hom_AMR"   "exac_non_TCGA.Hom_EAS"   "exac_non_TCGA.Hom_FIN"   "exac_non_TCGA.Hom_NFE"   "exac_non_TCGA.Hom_OTH"   "exac_non_TCGA.Hom_SAS"   "kgen.AC"                 "kgen.AF"                 "kgen.AFR_AF"             "kgen.AMR_AF"             "kgen.AN"                 "kgen.EAS_AF"             "kgen.EUR_AF"             "kgen.SAS_AF"             "SYMBOL"                  "Allele"                  "Existing_variation"      "Consequence"             "IMPACT"                  "CLIN_SIG"                "cDNA_position"           "CDS_position"            "Codons"                  "Protein_position"        "Amino_acids"             "DISTANCE"                "STRAND"                  "SYMBOL_SOURCE"           "Multiallelic"            "SIFT_call"               "SIFT_score"              "PolyPhen_call"           "PolyPhen_score"

# Check consistence of rownames and colnames

sum(rownames(gt.mx) != rownames(gq.mx))

## [1] 0

sum(rownames(gt.mx) != rownames(dp.mx))

## [1] 0

sum(rownames(gt.mx) != rownames(ref.mx))

## [1] 0

sum(rownames(gt.mx) != rownames(alt.mx))

## [1] 0

sum(colnames(gt.mx) != colnames(gq.mx))

## [1] 0

sum(colnames(gt.mx) != colnames(dp.mx))

## [1] 0

sum(colnames(gt.mx) != colnames(ref.mx))

## [1] 0

sum(colnames(gt.mx) != colnames(alt.mx))

## [1] 0

sum(rownames(gt.mx) != rownames(fixed.df))

## [1] 0

explore_data_before_filtering

Genotypes NA rates
Histogram of call rates per variant
Histograms of gq and dp in non-NA genotypes

gt.mx[1:3,c(1,2,541,542,738,739)]

##              100_S8_L007 101_S528_L008 9_S346_L008 HG00097 NA20828 NA20832
## Var000000001           0             0           0       0       0       0
## Var000000002           0             0           0       0       0       0
## Var000000003           0             0           0       0       0       0

# Prepare matrices with sub-groups data
gt_nfe.mx <- gt.mx[,542:739]
gq_nfe.mx <- gq.mx[,542:739]
dp_nfe.mx <- dp.mx[,542:739]
gt_wecare.mx <- gt.mx[,1:541]
gq_wecare.mx <- gq.mx[,1:541]
dp_wecare.mx <- dp.mx[,1:541]

gt_nfe.mx[1:5,1:5]

##              HG00097 HG00099 HG00100 HG00102 HG00106
## Var000000001       0       0       0       0       0
## Var000000002       0       0       0       0       0
## Var000000003       0       0       0       1       0
## Var000000004       1       2       1       1       2
## Var000000005       1       2       1       1       2

gt_wecare.mx[1:5,1:5]

##              100_S8_L007 101_S528_L008 102_S466_L008 103_S147_L007 104_S325_L008
## Var000000001           0             0             0             0             0
## Var000000002           0             0             0             0             0
## Var000000003           0             0             0             0             0
## Var000000004           0             2             2             2             1
## Var000000005           0             2             2             2             1

# --- Fractions of NA genotypes before filtering --- #

sum(is.na(gt.mx))/(nrow(gt.mx)*ncol(gt.mx)) # ~5%

## [1] 0.05399147

sum(is.na(gt_nfe.mx))/(nrow(gt_nfe.mx)*ncol(gt_nfe.mx)) # ~12%

## [1] 0.1254826

sum(is.na(gt_wecare.mx))/(nrow(gt_wecare.mx)*ncol(gt_wecare.mx)) # ~3%

## [1] 0.0278265

# --- Call rates per variant before filtering --- #

# All samples
x <- ncol(gt.mx)
y <- apply(gt.mx,1,function(z){1-sum(is.na(z))/x})
hist(y, breaks=50, main="Call rates per variant before filtering\n739 wecare-nfe samples", xlab="Call rates")

# nfe
x <- ncol(gt_nfe.mx)
y <- apply(gt_nfe.mx,1,function(z){1-sum(is.na(z))/x})
hist(y, breaks=50, main="Call rates per variant before filtering\n198 nfe samples", xlab="Call rates")

# wecare
x <- ncol(gt_wecare.mx)
y <- apply(gt_wecare.mx,1,function(z){1-sum(is.na(z))/x})
hist(y, breaks=50, main="Call rates per variant before filtering\n541 wecare samples", xlab="Call rates")

# --- gq  before filtering (when gt is not NA !) --- #

# All samples
hist(gq.mx[!is.na(gt.mx)], breaks=50, 
     main="Histogram of gq in non-NA genotypes before filtering\n739 wecare-nfe samples", xlab="gq")

hist(gq_nfe.mx[!is.na(gt_nfe.mx)], breaks=50, 
     main="Histogram of gq in non-NA genotypes before filtering\n198 nfe samples", xlab="gq")

hist(gq_wecare.mx[!is.na(gt_wecare.mx)], breaks=50, 
     main="Histogram of gq in non-NA genotypes before filtering\n541 wecare samples", xlab="gq")

# --- dp before filtering (when gt is not NA !) --- #

# All samples
hist(dp.mx[!is.na(gt.mx)], breaks=50, 
     main="Histogram of dp in non-NA genotypes before filtering\n739 wecare-nfe samples", xlab="dp")

hist(dp.mx[!is.na(gt.mx)], breaks=2500, xlim=c(0,100), 
     main="Histogram of dp in non-NA genotypes before filtering (zoom, 0:100)\n739 wecare-nfe samples", xlab="dp")

# nfe
hist(dp_nfe.mx[!is.na(gt_nfe.mx)], breaks=50, 
     main="Histogram of dp in non-NA genotypes before filtering\n198 nfe samples", xlab="dp")

hist(dp_nfe.mx[!is.na(gt_nfe.mx)], breaks=2500, xlim=c(0,100), 
     main="Histogram of dp in non-NA genotypes before filtering (zoom, 0:100)\n198 nfe samples", xlab="dp")

# wecare
hist(dp_wecare.mx[!is.na(gt_wecare.mx)], breaks=50, 
     main="Histogram of dp in non-NA genotypes before filtering\n541 wecare samples", xlab="dp")

hist(dp_wecare.mx[!is.na(gt_wecare.mx)], breaks=2500, xlim=c(0,100), 
     main="Histogram of dp in non-NA genotypes before filtering (zoom, 0:100)\n541 wecare samples", xlab="dp")

# Mean GQ and DP in non-NA genotypes before filtering
quantile(gq.mx[!is.na(gt.mx)], na.rm=TRUE) # median gq ~ 99

##   0%  25%  50%  75% 100% 
##    0   45   99   99   99

quantile(dp.mx[!is.na(gt.mx)], na.rm=TRUE) # median dp ~ 35

##   0%  25%  50%  75% 100% 
##    0   16   35   53 9901

mean(dp.mx[!is.na(gt.mx)], na.rm=TRUE) # mean dp ~89

## [1] 89.25717

# crude estimates for proportions of genotypes removed by filters
sum(dp.mx > max.dp, na.rm = T) / sum(!is.na(dp.mx))

## [1] 0.0122493

sum(dp.mx < min.dp, na.rm = T) / sum(!is.na(dp.mx))

## [1] 0.1797405

sum(gq.mx < min.gq, na.rm = T) / sum(!is.na(gq.mx))

## [1] 0.1190502

# Clean-up
rm(x,y, gt_nfe.mx, gq_nfe.mx, dp_nfe.mx, gt_wecare.mx, gq_wecare.mx, dp_wecare.mx)

filter_out_low_gq

Put NA to genotypes where gq < 20 : removes ~12% of non-NA genotypes

# num of genotypes to be removed (in !NA gt)
format(sum(gq.mx[!is.na(gt.mx)] < min.gq, na.rm=TRUE), big.mark=",") # ~687K

## [1] "686,697"

# Fraction of non-NA genotypes to be removed (in !NA gt)
sum(gq.mx[!is.na(gt.mx)] < min.gq, na.rm=TRUE)/sum(!is.na(gt.mx)) # ~12%

## [1] 0.1187019

# Apply filter (to gt only !)
NA -> gt.mx[ gq.mx < min.gq ]
NA -> gt_chr.mx[ gq.mx < min.gq ]

# Clean up
rm(min.gq)

explore_data_after_gq_filtering

# Prepare matrices with sub-groups data
gt_nfe.mx <- gt.mx[,542:739]
gq_nfe.mx <- gq.mx[,542:739]
dp_nfe.mx <- dp.mx[,542:739]
gt_wecare.mx <- gt.mx[,1:541]
gq_wecare.mx <- gq.mx[,1:541]
dp_wecare.mx <- dp.mx[,1:541]

# --- Fractions of NA genotypes after gq filtering --- #

sum(is.na(gt.mx))/(nrow(gt.mx)*ncol(gt.mx)) # ~17%

## [1] 0.1662845

sum(is.na(gt_nfe.mx))/(nrow(gt_nfe.mx)*ncol(gt_nfe.mx)) # ~22%

## [1] 0.225728

sum(is.na(gt_wecare.mx))/(nrow(gt_wecare.mx)*ncol(gt_wecare.mx)) # ~14%

## [1] 0.1445288

# --- Call rates per variant after gq filtering --- #

# All samples
x <- ncol(gt.mx)
y <- apply(gt.mx,1,function(z){1-sum(is.na(z))/x})
hist(y, breaks=50, main="Call rates per variant after gq filtering\n739 wecare-nfe samples", xlab="Call rates")

# nfe
x <- ncol(gt_nfe.mx)
y <- apply(gt_nfe.mx,1,function(z){1-sum(is.na(z))/x})
hist(y, breaks=50, main="Call rates per variant after gq filtering\n198 nfe samples", xlab="Call rates")

# wecare
x <- ncol(gt_wecare.mx)
y <- apply(gt_wecare.mx,1,function(z){1-sum(is.na(z))/x})
hist(y, breaks=50, main="Call rates per variant after gq filtering\n541 wecare samples", xlab="Call rates")

# --- gq  after gq filtering (when gt is not NA !) --- #

# All samples
hist(gq.mx[!is.na(gt.mx)], breaks=50, xlim=c(0,100), 
     main="Histogram of gq in non-NA genotypes after gq filtering\n739 wecare-nfe samples", xlab="gq")

hist(gq_nfe.mx[!is.na(gt_nfe.mx)], breaks=50, xlim=c(0,100), 
     main="Histogram of gq in non-NA genotypes after gq filtering\n198 nfe samples", xlab="gq")

hist(gq_wecare.mx[!is.na(gt_wecare.mx)], breaks=50, xlim=c(0,100), 
     main="Histogram of gq in non-NA genotypes after gq filtering\n541 wecare samples", xlab="gq")

# --- dp after gq filtering (when gt is not NA !) --- #

# All samples
hist(dp.mx[!is.na(gt.mx)], breaks=50, 
     main="Histogram of dp in non-NA genotypes after gq filtering\n739 wecare-nfe samples", xlab="dp")

hist(dp.mx[!is.na(gt.mx)], breaks=2500, xlim=c(0,100), 
     main="Histogram of dp in non-NA genotypes after gq filtering (zoom 0:100)\n739 wecare-nfe samples", xlab="dp")

# nfe
hist(dp_nfe.mx[!is.na(gt_nfe.mx)], breaks=50, 
     main="Histogram of dp in non-NA genotypes after gq filtering\n198 nfe samples", xlab="dp")

hist(dp_nfe.mx[!is.na(gt_nfe.mx)], breaks=2500, xlim=c(0,100), 
     main="Histogram of dp in non-NA genotypes after gq filtering (zoom 0:100)\n198 nfe samples", xlab="dp")

# wecare
hist(dp_wecare.mx[!is.na(gt_wecare.mx)], breaks=50, 
     main="Histogram of dp in non-NA genotypes after gq filtering\n541 wecare samples", xlab="dp")

hist(dp_wecare.mx[!is.na(gt_wecare.mx)], breaks=2500, xlim=c(0,100), 
     main="Histogram of dp in non-NA genotypes after gq filtering (zoom 0:100)\n541 wecare samples", xlab="dp")

# Mean GQ and DP in non-NA genotypes after gq filtering
mean(gq.mx[!is.na(gt.mx)], na.rm=TRUE) # ~ 81

## [1] 81.69107

mean(dp.mx[!is.na(gt.mx)], na.rm=TRUE) # ~ 100

## [1] 100.3238

# Clean-up
rm(x,y, gt_nfe.mx, gq_nfe.mx, dp_nfe.mx, gt_wecare.mx, gq_wecare.mx, dp_wecare.mx)

filter_out_low_dp

put NA to genotypes where dp < 10 : additionally removes ~5% of non-NA genotypes

# num of genotypes to be removed (in !NA gt after gq filtering)
format(sum(dp.mx[!is.na(gq.mx)] < min.dp, na.rm=TRUE), big.mark=",") # ~871K

## [1] "871,268"

# Fraction of genotypes to be removed (in !NA gt after gq filtering)
sum(dp.mx[!is.na(gq.mx)] < min.dp, na.rm=TRUE)/sum(!is.na(gt.mx)) # ~17%

## [1] 0.1708919

# Apply filter (to gt only !)
NA -> gt.mx[ dp.mx < min.dp ]
NA -> gt_chr.mx[ dp.mx < min.dp ]

# Clean up
rm(min.dp)

explore_data_after_gq_min-dp_filtering

# Prepare matrices with sub-groups data
gt_nfe.mx <- gt.mx[,542:739]
gq_nfe.mx <- gq.mx[,542:739]
dp_nfe.mx <- dp.mx[,542:739]
gt_wecare.mx <- gt.mx[,1:541]
gq_wecare.mx <- gq.mx[,1:541]
dp_wecare.mx <- dp.mx[,1:541]

# --- Fractions of NA genotypes after gq filtering --- #

sum(is.na(gt.mx))/(nrow(gt.mx)*ncol(gt.mx)) # ~21%

## [1] 0.2075736

sum(is.na(gt_nfe.mx))/(nrow(gt_nfe.mx)*ncol(gt_nfe.mx)) # ~25%

## [1] 0.2469261

sum(is.na(gt_wecare.mx))/(nrow(gt_wecare.mx)*ncol(gt_wecare.mx)) # ~19%

## [1] 0.193171

# --- Call rates per variant after gq-dp filtering --- #

# All samples
x <- ncol(gt.mx)
y <- apply(gt.mx,1,function(z){1-sum(is.na(z))/x})
hist(y, breaks=50, main="Call rates per variant after gq-low-dp filtering\n739 wecare-nfe samples", xlab="Call rates")

# nfe
x <- ncol(gt_nfe.mx)
y <- apply(gt_nfe.mx,1,function(z){1-sum(is.na(z))/x})
hist(y, breaks=50, main="Call rates per variant after gq-low-dp filtering\n198 nfe samples", xlab="Call rates")

# wecare
x <- ncol(gt_wecare.mx)
y <- apply(gt_wecare.mx,1,function(z){1-sum(is.na(z))/x})
hist(y, breaks=50, main="Call rates per variant after gq-low-dp filtering\n541 wecare samples", xlab="Call rates")

# --- gq  after gq-dp filtering (when gt is not NA !) --- #

# All samples
hist(gq.mx[!is.na(gt.mx)], breaks=50, xlim=c(0,100), 
     main="Histogram of gq in non-NA genotypes after gq-low-dp filtering\n739 wecare-nfe samples", xlab="gq")

hist(gq_nfe.mx[!is.na(gt_nfe.mx)], breaks=50, xlim=c(0,100), 
     main="Histogram of gq in non-NA genotypes after gq-low-dp filtering\n198 nfe samples", xlab="gq")

hist(gq_wecare.mx[!is.na(gt_wecare.mx)], breaks=50, xlim=c(0,100), 
     main="Histogram of gq in non-NA genotypes after gq-low-dp filtering\n541 wecare samples", xlab="gq")

# --- dp after gq-dp filtering (when gt is not NA !) --- #

# All samples
hist(dp.mx[!is.na(gt.mx)], breaks=50, 
     main="Histogram of dp in non-NA genotypes after gq-low-dp filtering\n739 wecare-nfe samples", xlab="dp")

hist(dp.mx[!is.na(gt.mx)], breaks=2500, xlim=c(0,100), 
     main="Histogram of dp in non-NA genotypes after gq-low-dp filtering (zoom, 0:100)\n739 wecare-nfe samples", xlab="dp")

# nfe
hist(dp_nfe.mx[!is.na(gt_nfe.mx)], breaks=50, 
     main="Histogram of dp in non-NA genotypes after gq-low-dp filtering\n198 nfe samples", xlab="dp")

hist(dp_nfe.mx[!is.na(gt_nfe.mx)], breaks=250, xlim=c(0,100), 
     main="Histogram of dp in non-NA genotypes after gq-low-dp filtering (zoom, 0:100)\n198 nfe samples", xlab="dp")

# wecare
hist(dp_wecare.mx[!is.na(gt_wecare.mx)], breaks=50, 
     main="Histogram of dp in non-NA genotypes after gq-low-dp filtering\n541 wecare samples", xlab="dp")

hist(dp_wecare.mx[!is.na(gt_wecare.mx)], breaks=2500, xlim=c(0,100), 
     main="Histogram of dp in non-NA genotypes after gq-low-dp filtering (zoom, 0:100)\n541 wecare samples", xlab="dp")

# Mean GQ and DP in non-NA genotypes after gq-dp filtering
mean(gq.mx[!is.na(gt.mx)], na.rm=TRUE) # ~ 85

## [1] 84.50231

mean(dp.mx[!is.na(gt.mx)], na.rm=TRUE) # ~ 105

## [1] 105.1593

# Clean-up
rm(x,y, gt_nfe.mx, gq_nfe.mx, dp_nfe.mx, gt_wecare.mx, gq_wecare.mx, dp_wecare.mx)

filter_out_high_dp

put NA to genotypes where dp > 1000 : removes <<1% of non-NA genotypes

# num of genotypes to be removed (in !NA gt after previous filtering)
format(sum(dp.mx[!is.na(gq.mx)] > max.dp, na.rm=TRUE), big.mark=",") # ~73K

## [1] "72,630"

# Fraction of genotypes to be removed (in !NA gt after previous filtering)
sum(dp.mx[!is.na(gq.mx)] > max.dp, na.rm=TRUE)/sum(!is.na(gt.mx)) # ~1%

## [1] 0.01498803

# Apply filter (to gt only !)
NA -> gt.mx[ dp.mx > max.dp ]
NA -> gt_chr.mx[ dp.mx > max.dp ]

# Clean up
rm(max.dp)

explore_data_after_gq_dp_filtering

# Prepare matrices with sub-groups data
gt_nfe.mx <- gt.mx[,542:739]
gq_nfe.mx <- gq.mx[,542:739]
dp_nfe.mx <- dp.mx[,542:739]
gt_wecare.mx <- gt.mx[,1:541]
gq_wecare.mx <- gq.mx[,1:541]
dp_wecare.mx <- dp.mx[,1:541]

# --- Fractions of NA genotypes after gq-dp filtering --- #

sum(is.na(gt.mx))/(nrow(gt.mx)*ncol(gt.mx)) # ~22%

## [1] 0.2194281

sum(is.na(gt_nfe.mx))/(nrow(gt_nfe.mx)*ncol(gt_nfe.mx)) # ~25%

## [1] 0.2469261

sum(is.na(gt_wecare.mx))/(nrow(gt_wecare.mx)*ncol(gt_wecare.mx)) # ~20%

## [1] 0.2093641

# --- Call rates per variant after gq-dp filtering --- #

# All samples
x <- ncol(gt.mx)
y <- apply(gt.mx,1,function(z){1-sum(is.na(z))/x})
hist(y, breaks=50, main="Call rates per variant after gq-dp filtering\n739 wecare-nfe samples", xlab="Call rates")

# nfe
x <- ncol(gt_nfe.mx)
y <- apply(gt_nfe.mx,1,function(z){1-sum(is.na(z))/x})
hist(y, breaks=50, main="Call rates per variant after gq-dp filtering\n198 nfe samples", xlab="Call rates")

# wecare
x <- ncol(gt_wecare.mx)
y <- apply(gt_wecare.mx,1,function(z){1-sum(is.na(z))/x})
hist(y, breaks=50, main="Call rates per variant after gq-dp filtering\n541 wecare samples", xlab="Call rates")

# --- gq  after gq-dp filtering (when gt is not NA !) --- #

# All samples
hist(gq.mx[!is.na(gt.mx)], breaks=50, xlim=c(0,100), 
     main="Histogram of gq in non-NA genotypes after gq-dp filtering\n739 wecare-nfe samples", xlab="gq")

hist(gq_nfe.mx[!is.na(gt_nfe.mx)], breaks=50, xlim=c(0,100), 
     main="Histogram of gq in non-NA genotypes after gq-dp filtering\n198 nfe samples", xlab="gq")

hist(gq_wecare.mx[!is.na(gt_wecare.mx)], breaks=50, xlim=c(0,100), 
     main="Histogram of gq in non-NA genotypes after gq-dp filtering\n541 wecare samples", xlab="gq")

# --- dp after gq-dp filtering (when gt is not NA !) --- #

# All samples
hist(dp.mx[!is.na(gt.mx)], breaks=50, 
     main="Histogram of dp in non-NA genotypes after gq-dp filtering\n739 wecare-nfe samples", xlab="dp")

hist(dp.mx[!is.na(gt.mx)], breaks=250, xlim=c(0,100), 
     main="Histogram of dp in non-NA genotypes after gq-dp filtering (zoom, 0:100)\n739 wecare-nfe samples", xlab="dp")

# nfe
hist(dp_nfe.mx[!is.na(gt_nfe.mx)], breaks=50, 
     main="Histogram of dp in non-NA genotypes after gq-dp filtering\n198 nfe samples", xlab="dp")

hist(dp_nfe.mx[!is.na(gt_nfe.mx)], breaks=250, xlim=c(0,100), 
     main="Histogram of dp in non-NA genotypes after gq-dp filtering (zoom, 0:100)\n198 nfe samples", xlab="dp")

# wecare
hist(dp_wecare.mx[!is.na(gt_wecare.mx)], breaks=50, 
     main="Histogram of dp in non-NA genotypes after gq-dp filtering\n541 wecare samples", xlab="dp")

hist(dp_wecare.mx[!is.na(gt_wecare.mx)], breaks=250, xlim=c(0,100), 
     main="Histogram of dp in non-NA genotypes after gq-dp filtering (zoom, 0:100)\n541 wecare samples", xlab="dp")

# Mean GQ and DP in non-NA genotypes after gq-dp filtering
mean(gq.mx[!is.na(gt.mx)], na.rm=TRUE) # ~ 84

## [1] 84.28232

mean(dp.mx[!is.na(gt.mx)], na.rm=TRUE) # ~ 73

## [1] 73.08944

# Clean-up
rm(x,y, gt_nfe.mx, gq_nfe.mx, dp_nfe.mx, gt_wecare.mx, gq_wecare.mx, dp_wecare.mx)

Filtering by fraction of reads supporting ALt allele

Check Alt allele fraction in Hom-Ref, Het and Hom-Alt

# Prepare alt_fraction.mx
alt_fraction.mx <- alt.mx / (alt.mx + ref.mx)
dim(alt_fraction.mx)

## [1] 8275  739

sum(is.na(alt_fraction.mx))

## [1] 336781

sum(is.na(gt.mx))

## [1] 1341852

# Only consider alt_fraction, where genotypes are not NA
NA -> alt_fraction.mx[is.na(gt.mx)]

sum(is.na(alt_fraction.mx))

## [1] 1348217

sum(is.na(gt.mx))

## [1] 1341852

# Explore alt_fraction in hom.ref

alt_fraction_hom_ref.mx <- alt_fraction.mx
NA -> alt_fraction_hom_ref.mx[gt.mx != 0]
sum(!is.na(alt_fraction_hom_ref.mx))

## [1] 4530694

hist(alt_fraction_hom_ref.mx, xlim=c(0,1))

sum(alt_fraction_hom_ref.mx > 0.2, na.rm = T)

## [1] 81

# Explore alt_fraction in het

alt_fraction_het.mx <- alt_fraction.mx
NA -> alt_fraction_het.mx[gt.mx != 1]
sum(!is.na(alt_fraction_het.mx))

## [1] 154851

hist(alt_fraction_het.mx, xlim=c(0,1))

sum(alt_fraction_het.mx < 0.2, na.rm = T)

## [1] 5988

sum(alt_fraction_het.mx > 0.8, na.rm = T)

## [1] 1341

# Explore alt_fraction in hom.alt

alt_fraction_hom_alt.mx <- alt_fraction.mx
NA -> alt_fraction_hom_alt.mx[gt.mx != 2]
sum(!is.na(alt_fraction_hom_alt.mx))

## [1] 81463

hist(alt_fraction_hom_alt.mx, xlim=c(0,1))

sum(alt_fraction_hom_alt.mx < 0.2, na.rm = T)

## [1] 31

# Clean-up (keep matrices needed in teh next chunk!)
rm(alt_fraction.mx)

Update genotypes matrix

NA -> gt.mx[alt_fraction_hom_ref.mx > 0.2]
NA -> gt_chr.mx[alt_fraction_hom_ref.mx > 0.2]

NA -> gt.mx[alt_fraction_het.mx > 0.8]
NA -> gt_chr.mx[alt_fraction_het.mx > 0.8]

NA -> gt.mx[alt_fraction_het.mx < 0.2]
NA -> gt_chr.mx[alt_fraction_het.mx < 0.2]

NA -> gt.mx[alt_fraction_hom_alt.mx < 0.8]
NA -> gt_chr.mx[alt_fraction_hom_alt.mx < 0.8]

rm(alt_fraction_hom_ref.mx, alt_fraction_het.mx, alt_fraction_hom_alt.mx)

Repeat assessment of allelic fractions in Hom-Ref, Het and Hom-Alt

# Prepare alt_fraction.mx
alt_fraction.mx <- alt.mx / (alt.mx + ref.mx)
dim(alt_fraction.mx)

## [1] 8275  739

sum(is.na(alt_fraction.mx))

## [1] 336781

sum(is.na(gt.mx))

## [1] 1349295

# Only consider alt_fraction, where genotypes are not NA
NA -> alt_fraction.mx[is.na(gt.mx)]

sum(is.na(alt_fraction.mx))

## [1] 1355660

sum(is.na(gt.mx))

## [1] 1349295

# Explore alt_fraction in hom.ref

alt_fraction_hom_ref.mx <- alt_fraction.mx
NA -> alt_fraction_hom_ref.mx[gt.mx != 0]
sum(!is.na(alt_fraction_hom_ref.mx))

## [1] 4530613

hist(alt_fraction_hom_ref.mx, xlim=c(0,1))

sum(alt_fraction_hom_ref.mx > 0.2, na.rm = T)

## [1] 0

# Explore alt_fraction in het

alt_fraction_het.mx <- alt_fraction.mx
NA -> alt_fraction_het.mx[gt.mx != 1]
sum(!is.na(alt_fraction_het.mx))

## [1] 147522

hist(alt_fraction_het.mx, xlim=c(0,1))

sum(alt_fraction_het.mx < 0.2, na.rm = T)

## [1] 0

sum(alt_fraction_het.mx > 0.8, na.rm = T)

## [1] 0

# Explore alt_fraction in hom.alt

alt_fraction_hom_alt.mx <- alt_fraction.mx
NA -> alt_fraction_hom_alt.mx[gt.mx != 2]
sum(!is.na(alt_fraction_hom_alt.mx))

## [1] 81430

hist(alt_fraction_hom_alt.mx, xlim=c(0,1))

sum(alt_fraction_hom_alt.mx < 0.2, na.rm = T)

## [1] 0

# Clean-up
rm(alt_fraction.mx, alt_fraction_hom_ref.mx, alt_fraction_het.mx, 
   alt_fraction_hom_alt.mx, ref.mx, alt.mx)

explore_data_after_gq_dp_af_filtering

# Prepare matrices with sub-groups data
gt_nfe.mx <- gt.mx[,542:739]
gq_nfe.mx <- gq.mx[,542:739]
dp_nfe.mx <- dp.mx[,542:739]
gt_wecare.mx <- gt.mx[,1:541]
gq_wecare.mx <- gq.mx[,1:541]
dp_wecare.mx <- dp.mx[,1:541]

# --- Fractions of NA genotypes after gq-dp-af filtering --- #

sum(is.na(gt.mx))/(nrow(gt.mx)*ncol(gt.mx)) # ~22%

## [1] 0.2206452

sum(is.na(gt_nfe.mx))/(nrow(gt_nfe.mx)*ncol(gt_nfe.mx)) # ~25%

## [1] 0.2472721

sum(is.na(gt_wecare.mx))/(nrow(gt_wecare.mx)*ncol(gt_wecare.mx)) # ~21%

## [1] 0.2109

# --- Call rates per variant after gq-dp-af filtering --- #

# All samples
x <- ncol(gt.mx)
y <- apply(gt.mx,1,function(z){1-sum(is.na(z))/x})
hist(y, breaks=50, main="Call rates per variant after gq-dp-af filtering\n739 wecare-nfe samples", xlab="Call rates")

# nfe
x <- ncol(gt_nfe.mx)
y <- apply(gt_nfe.mx,1,function(z){1-sum(is.na(z))/x})
hist(y, breaks=50, main="Call rates per variant after gq-dp-af filtering\n198 nfe samples", xlab="Call rates")

# wecare
x <- ncol(gt_wecare.mx)
y <- apply(gt_wecare.mx,1,function(z){1-sum(is.na(z))/x})
hist(y, breaks=50, main="Call rates per variant after gq-dp-af filtering\n541 wecare samples", xlab="Call rates")

# --- gq  after gq-dp-af filtering (when gt is not NA !) --- #

# All samples
hist(gq.mx[!is.na(gt.mx)], breaks=50, xlim=c(0,100), 
     main="Histogram of gq in non-NA genotypes after gq-dp-af filtering\n739 wecare-nfe samples", xlab="gq")

hist(gq_nfe.mx[!is.na(gt_nfe.mx)], breaks=50, xlim=c(0,100), 
     main="Histogram of gq in non-NA genotypes after gq-dp-af filtering\n198 nfe samples", xlab="gq")

hist(gq_wecare.mx[!is.na(gt_wecare.mx)], breaks=50, xlim=c(0,100), 
     main="Histogram of gq in non-NA genotypes after gq-dp-af filtering\n541 wecare samples", xlab="gq")

# --- dp after gq-dp-af filtering (when gt is not NA !) --- #

# All samples
hist(dp.mx[!is.na(gt.mx)], breaks=50, 
     main="Histogram of dp in non-NA genotypes after gq-dp-af filtering\n739 wecare-nfe samples", xlab="dp")

hist(dp.mx[!is.na(gt.mx)], breaks=250, xlim=c(0,100), 
     main="Histogram of dp in non-NA genotypes after gq-dp-af filtering (zoom, 0:100)\n739 wecare-nfe samples", xlab="dp")

# nfe
hist(dp_nfe.mx[!is.na(gt_nfe.mx)], breaks=50, 
     main="Histogram of dp in non-NA genotypes after gq-dp-af filtering\n198 nfe samples", xlab="dp")

hist(dp_nfe.mx[!is.na(gt_nfe.mx)], breaks=250, xlim=c(0,100), 
     main="Histogram of dp in non-NA genotypes after gq-dp-af filtering (zoom, 0:100)\n198 nfe samples", xlab="dp")

# wecare
hist(dp_wecare.mx[!is.na(gt_wecare.mx)], breaks=50, 
     main="Histogram of dp in non-NA genotypes after gq-dp-af filtering\n541 wecare samples", xlab="dp")

hist(dp_wecare.mx[!is.na(gt_wecare.mx)], breaks=250, xlim=c(0,100), 
     main="Histogram of dp in non-NA genotypes after gq-dp-af filtering (zoom, 0:100)\n541 wecare samples", xlab="dp")

# Mean GQ and DP in non-NA genotypes after gq-dp-af filtering
mean(gq.mx[!is.na(gt.mx)], na.rm=TRUE) # ~ 84

## [1] 84.31422

mean(dp.mx[!is.na(gt.mx)], na.rm=TRUE) # ~ 73

## [1] 73.09395

# Clean-up
rm(x,y, gt_nfe.mx, gq_nfe.mx, dp_nfe.mx, gt_wecare.mx, gq_wecare.mx, dp_wecare.mx)

filter_variants_by_final_call_rate

Remove variants with call rate < 85% in either nfe or wecare after the genotypes filtering:
removes 6,065 (~73%) variants (8,275 -> 2,210)
Most of the variants were removed because of low call rate in NFE

dim(gt.mx)

## [1] 8275  739

# Estimate callrates in WECARE samples
wecare_cases <- c(rep(T,541),rep(F,198))
length(wecare_cases)

## [1] 739

sum(wecare_cases)

## [1] 541

x_wecare <- ncol(gt.mx[,wecare_cases])
y_wecare <- apply(gt.mx[,wecare_cases], 1, function(z){1-sum(is.na(z))/x_wecare})
length(y_wecare)

## [1] 8275

y_wecare[1:7]

## Var000000001 Var000000002 Var000000003 Var000000004 Var000000005 Var000000006 Var000000010 
##    0.9889094    0.9889094    0.9426987    0.7338262    0.5730129    0.7707948    0.9020333

var_wecare.retained <- y_wecare >= min.call.rate
sum(var_wecare.retained) # 4,916 to be retained

## [1] 4916

1 - sum(var_wecare.retained)/nrow(gt.mx[,wecare_cases]) # ~40% to be removed

## [1] 0.4059215

# Estimate callrates in NFE samples
nfe_cases <- !wecare_cases
sum(nfe_cases)

## [1] 198

x_nfe <- ncol(gt.mx[,nfe_cases])
y_nfe <- apply(gt.mx[,nfe_cases], 1, function(z){1-sum(is.na(z))/x_nfe})
length(y_nfe)

## [1] 8275

y_nfe[1:7]

## Var000000001 Var000000002 Var000000003 Var000000004 Var000000005 Var000000006 Var000000010 
##    0.3636364    0.8282828    0.8585859    0.9090909    0.9090909    0.9141414    0.5151515

var_nfe.retained <- y_nfe >= min.call.rate
sum(var_nfe.retained) # 3,565 to be retained

## [1] 3565

1 - sum(var_nfe.retained)/nrow(gt.mx[,nfe_cases]) # ~57% to be removed

## [1] 0.5691843

# Estimate the proportion of variants to be removed
var.retained <- var_wecare.retained & var_nfe.retained
sum(var.retained) # 2,210 to be retained

## [1] 2210

1 - sum(var.retained)/nrow(gt.mx) # ~73% to be removed

## [1] 0.7329305

# Remove variants with loaw call rates
gt.mx <- gt.mx[ var.retained, ]
gt_chr.mx <- gt_chr.mx[ var.retained, ]

dp.mx <- dp.mx[ var.retained, ]
gq.mx <- gq.mx[ var.retained, ]
fixed.df <- fixed.df[ var.retained, ]

# Clean-up
rm(min.call.rate, var.retained, x_wecare, y_wecare, x_nfe, y_nfe, 
   wecare_cases, nfe_cases, var_wecare.retained, var_nfe.retained)

remove_variants_with_the_uniform_genotypes_accross_all_samples

Remove 324 variants: 2,210 -> 1,886
In some variants the diverse genotypes were removed during the above filtering

# Check that there is no all-NA variants
non_NA_count.udf <- function(x){sum(!is.na(x))}
all_NA <- apply(gt.mx, 1, non_NA_count.udf) == 0
sum(all_NA) # 0

## [1] 0

# Function to detect uniform numeric vector
uniform_vector.udf <- function(x){
  if(min(x, na.rm=TRUE) == max(x, na.rm=TRUE)){return(TRUE)} else {return(FALSE)}}

# Variants with uniform genotypes accross all samples 
uniform_genotypes <- apply(gt.mx, 1, uniform_vector.udf)
summary(uniform_genotypes)

##    Mode   FALSE    TRUE 
## logical    1886     324

sum(uniform_genotypes) # 324

## [1] 324

# Remove variants with uniform genotypes accross all samples
gt.mx <- gt.mx[!uniform_genotypes,]
gt_chr.mx <- gt_chr.mx[!uniform_genotypes,]
gq.mx <- gq.mx[!uniform_genotypes,]
dp.mx <- dp.mx[!uniform_genotypes,]

fixed.df <- fixed.df[!uniform_genotypes,]

dim(gt.mx)

## [1] 1886  739

dim(gt_chr.mx)

## [1] 1886  739

dim(gq.mx)

## [1] 1886  739

dim(dp.mx)

## [1] 1886  739

dim(fixed.df)

## [1] 1886  110

# Clean-up
rm(non_NA_count.udf, all_NA, uniform_vector.udf, uniform_genotypes)

explore_data_after_filtering